

initialBulkComps = [0.894842208	0.059505645	0.118
0.905480461	0.060773481	0.05
0.900368082	0.067929065	0.161
0.894445082	0.071090047	0.124
0.889997816	0.078534031	0.183
0.904827309	0.079896907	0.11
0.892348295	0.082446809	0.16
0.894052698	0.08275463	0.13
0.899477263	0.082840237	0.174
0.891383871	0.088685015	0.11
0.889997816	0.090002493	0.21
0.892001701	0.093333333	0.15
0.897215904	0.093499093	0.217
0.874501581	0.093924943	0.33
0.895722335	0.09439528	0.11
0.893265325	0.098756029	0.201
0.909774831	0.099437148	0.1
0.899298855	0.10184667	0.181
0.899368602	0.101872421	0.166
0.888491133	0.103641457	0.15
0.899333851	0.103758432	0.158
0.900981012	0.105968972	0.17
0.88676188	0.110974862	0.16
0.896143917	0.119318182	0.2
0.868410258	0.126005362	0.25
0.898557508	0.139534884	0.19
0.919081956	0.140625	0.07
0.892001701	0.171446154	0.1485
0.892001701	0.171446154	0.1485
0.905480461	0.214249031	0.0495
0.930783185	0.246376812	0.054]; 
%%
figure(1)
close
figure(1)
hold on
FigureTitle = 'FigureS2-Bulk-MeltConditions';
set(gcf,'name',regexprep(FigureTitle,'\_*',''))

axesMajors = [0 10 0 30; 0 10 0 30 ;...
    0 10 0 30 ];




numcolumns = 3;
numrows = 3;


things2plot_labels = {'Mg#(mol) Mantle Source' 'NaK#(wt%) Mantle Source' 'TiO_2(wt%) Mantle Source'};


for kk = 1:3
        figure(1)
        set(gcf, 'Units', 'Inches', 'Position', [8.0556 3.8889 10.7083 9.7778], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        %set(gcf, 'Units', 'Inches', 'Position',  [3.7917 7 16.4861 5.6111], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        subaxis(numrows,numcolumns,kk,'Spacing',.03,'Margin',0.015,'MarginLeft',.05,'Padding',0.00,'PaddingLeft',0.015,'PaddingBottom',.05)
        
       % axis(axesMajors(kk,:))
        hold on
        axis square
        
       
        for iik=1:numberModels

            
            layerplot_color_here = eval(char(layerplot_color(iik)));
            layerplot_color_here_fill = eval(char(layerplot_color_fill(iik)));
            layerplot_color_here_edge = eval(char(layerplot_color_edge(iik)));
             mymarkersize = layerplot_markersize(iik);
             linewidth_here = layerplot_linewidth(iik);

          
           plot(initialBulkComps(:,kk),...
            eval(sprintf('%s_GSBoundary', TabNames{iik})), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',mymarkersize,'LineWidth',linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge)


            Highlight = [8 18];
            plot(initialBulkComps(Highlight,kk),...
            eval(sprintf('%s_GSBoundary(%s,:)', TabNames{iik},'Highlight')), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here,...
            'MarkerSize',1.5*mymarkersize,'LineWidth',1.5*linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge,'handlevisibility','off')
       

   
            %text(2,25,sprintf('+%g ppm',WaterTarget(kk)),'FontName','Times New Roman','FontSize',17)
        
        end
        
        ticklengthgood = 4*[0.01 0.025]; 
        set(gca,'ticklength',ticklengthgood)
        set(gca,'fontsize', 13,'LineWidth',0.7,'FontName','Times New Roman')
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
      
       
        ax=gca;

         
         xlabel(things2plot_labels{kk})
         
         if kk==1
         ylabel('Gar-Sp Boundary (kbar)')
          legend([names4legend],'FontSize',10)
         end
        
        
        box on
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
     
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
end



    for kk = [1:3]

        subaxis(numrows,numcolumns,kk+3)
        
       % axis(axesMajors(kk,:))
        hold on
        axis square
        
       
        for iik=1:numberModels
            
            
           layerplot_color_here = eval(char(layerplot_color(iik)));
            layerplot_color_here_fill = eval(char(layerplot_color_fill(iik)));
            layerplot_color_here_edge = eval(char(layerplot_color_edge(iik)));
             mymarkersize = layerplot_markersize(iik);
             linewidth_here = layerplot_linewidth(iik);
  
            


           plot(initialBulkComps(:,kk),...
            eval(sprintf('%s_MeltExtractP', TabNames{iik})), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',mymarkersize,'LineWidth',linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge)
   
        
             plot(initialBulkComps(:,kk),...
            eval(sprintf('%s_LastMeltP', TabNames{iik})), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',mymarkersize,'LineWidth',linewidth_here,'MarkerEdgeColor',fuschia)
        
        
        
            plot(initialBulkComps(:,kk),...
            eval(sprintf('%s_GSBoundary', TabNames{iik})), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',5,'LineWidth',layerplot_linewidth(iik),'MarkerEdgeColor','k')


            plot(initialBulkComps(Highlight,kk),...
            eval(sprintf('%s_MeltExtractP(%s,:)', TabNames{iik},'Highlight')), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',1.5*mymarkersize,'LineWidth',2*linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge,'HandleVisibility','off')
   
          
        end
        
        ticklengthgood = 4*[0.01 0.025]; 
        set(gca,'ticklength',ticklengthgood)
        set(gca,'fontsize', 13,'LineWidth',0.7,'FontName','Times New Roman')
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
         
       
        ax=gca;

         
         xlabel(things2plot_labels{kk})
         
         if kk==1
         ylabel('Extrac. and G-S Pressure (kbar)')
         end
        
        
        box on
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
     
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
    end


    
    for kk = [1:3]
       subaxis(numrows,numcolumns,kk+6)
        
       % axis(axesMajors(kk,:))
        hold on
        axis square
        
       
        for iik=1:numberModels
      
            
                 
           layerplot_color_here = eval(char(layerplot_color(iik)));
            layerplot_color_here_fill = eval(char(layerplot_color_fill(iik)));
            layerplot_color_here_edge = eval(char(layerplot_color_edge(iik)));
             mymarkersize = layerplot_markersize(iik);
             linewidth_here = layerplot_linewidth(iik);
            
          

           plot(initialBulkComps(:,kk),...
            eval(sprintf('%s_Tdiff', TabNames{iik})), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',layerplot_markersize(iik),'LineWidth',linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge)


            %text(2,25,sprintf('+%g ppm',WaterTarget(kk)),'FontName','Times New Roman','FontSize',17)
        
            plot(initialBulkComps(Highlight,kk),...
            eval(sprintf('%s_Tdiff(%s,:)', TabNames{iik},'Highlight')), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',1.5*layerplot_markersize(iik),'LineWidth',1.5*linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge,'HandleVisibility','off')
          
        
          
        end
        
        ticklengthgood = 4*[0.01 0.025]; 
        set(gca,'ticklength',ticklengthgood)
        set(gca,'fontsize', 13,'LineWidth',0.7,'FontName','Times New Roman')
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        
       
        ax=gca;

         
         xlabel(things2plot_labels{kk})
         
         if kk==1
         ylabel('Gar-Sp \DeltaT(\circC)')
         end
        
        
        box on
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
     
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
    end

    
      
        
axes=[.86 .94 21 28
    0.05 .25 21 28
    0 .4 21 28];
for ll = 1:3
subaxis(3,3,ll)
axis(axes(ll,:))
box on
axis square


 ax = gca;
 ax.XAxis.MinorTick = 'on';
increment =  findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
if ll==3 
    increment=.05;
end
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment=.5; 
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);

set(gca,'YDir','Reverse')
end


 axes=[.86 .94 10 35
    0.05 .25 10 35
    0 .4 10 35];


for ll = 1:3
subaxis(3,3,ll+3)
axis(axes(ll,:))
end

for ll = 4:6
subaxis(3,3,ll)
box on
axis square
 ax = gca;
 ax.XAxis.MinorTick = 'on';
increment =  findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
if ll==6
    increment=.05;
end
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
increment=1; 
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
set(gca,'YDir','Reverse')
end



for ll = 6:9
subaxis(3,3,ll)

axis square
% set(gca,'fontsize', 15,'LineWidth',1)
% set(gca,'XColor', 'k')
% set(gca,'YColor', 'k')

 ax = gca;

 ax.XAxis.MinorTick = 'on';
increment =  findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
if ll==9 
    increment=.05;
end
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment =5; % findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
box on
set(gca,'YDir','Normal')
end


%%

figure(2)
close
figure(2)
hold on
FigureTitle = 'FigureS3B-Bulk-Zeta';
set(gcf,'name',regexprep(FigureTitle,'\_*',''))

axesMajors = [0 10 0 30; 0 10 0 30 ;...
    0 10 0 30 ];




numcolumns = 3;
numrows = 1;


things2plot_labels = {'Mg#(mol) Mantle Source' 'NaK#(wt%) Mantle Source' 'TiO_2(wt%) Mantle Source'};

for kk = [1:3]


 
        %for 3x2
        %set(gcf, 'Units', 'Inches', 'Position', [8.9306 0 6.8333 9.7778], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        % set(gcf, 'Units', 'Inches', 'Position',[6.5278 0.1667 8.8056 12.7222], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        % subaxis(numrows,numcolumns,n,'Spacing',.08,'Margin',0.02,'MarginBottom',.06,'MarginLeft',.06,'Padding',0,'PaddingLeft',0.02)
        
        %for 3x3
        %set(gcf, 'Units', 'Inches', 'Position', [8.0556 3.8889 10.7083 9.7778], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        set(gcf, 'Units', 'Inches', 'Position',  [3.7917 7 16.4861 5.6111], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        subaxis(numrows,numcolumns,kk,'Spacing',.03,'Margin',0.015,'MarginLeft',.05,'Padding',0.00,'PaddingLeft',0.015,'PaddingBottom',.05)
        
     
        
       % axis(axesMajors(kk,:))
        hold on
        axis square
        
       
        for iik=1:numberModels

            
            layerplot_color_here = eval(char(layerplot_color(iik)));
            layerplot_color_here_fill = eval(char(layerplot_color_fill(iik)));
            layerplot_color_here_edge = eval(char(layerplot_color_edge(iik)));
             mymarkersize = layerplot_markersize(iik);
             linewidth_here = layerplot_linewidth(iik);
            
          
            
            
          

            plot(initialBulkComps(:,kk),...
            eval(sprintf('%s_WGP', TabNames{iik})), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',layerplot_markersize(iik),'LineWidth',linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge)
        
        
            Highlight = [8 18]; 
            plot(initialBulkComps(Highlight,kk),...
            eval(sprintf('%s_WGP(%s,:)', TabNames{iik},'Highlight')), ...
            char(layerplot_marker(iik)), 'Color',layerplot_color_here, 'MarkerFaceColor',layerplot_color_here_fill,...
            'MarkerSize',1.5*layerplot_markersize(iik),'LineWidth',1.5.*linewidth_here,'MarkerEdgeColor',layerplot_color_here_edge,'HandleVisibility','off')


            %text(2,25,sprintf('+%g ppm',WaterTarget(kk)),'FontName','Times New Roman','FontSize',17)
        
        
          
        
          
        end
        
        ticklengthgood = 4*[0.01 0.025]; 
        set(gca,'ticklength',ticklengthgood)
        set(gca,'fontsize', 18,'LineWidth',0.7,'FontName','Times New Roman')
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
  
       
        ax=gca;
       ytickformat(ax, 'percentage');
         
         xlabel(things2plot_labels{kk})
         
         if kk==1
         ylabel('\zeta_P')
          legend(names4legend,'FontSize',10)
         end
        
        
        box on
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
     
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
        
   
     
end



axes=[.86 .94 0 50
    0.05 .25 0 50
    0 .4 0 50];
axis square

for data2plot_num = 1:3
subaxis(1,3,data2plot_num)
axis(axes(data2plot_num,:))
box on
axis square


ax = gca;

ax.XAxis.MinorTick = 'on';
increment =  findBestIncrement(ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
if data2plot_num ==6
    increment = .05;
end

if data2plot_num ==3
    increment = .05;
end
ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
ax.YAxis.MinorTick = 'on';
increment = findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);

end